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Abstract 

By extending a dynamical mean-field approximation (DMA) previously pro- 
posed by the author [H. Hasegawa, Phys. Rev. E 67, 41903 (2003)], we have 
developed a semianalytical theory which takes into account a wide range of 
couplings in a small-world network. Our network consists of noisy A^-unit 
FitzHugh-Nagumo (FN) neurons with couplings whose average coordination 
number Z may change from local (Z <C A^) to global couplings (Z = N — 1) 
and/or whose concentration of random couplings p is allowed to vary from 
regular (p = 0) to completely random (p=l). We have taken into account 
three kinds of spatial correlations: the on-site correlation, the correlation 
for a coupled pair and that for a pair without direct couplings. The orig- 
inal 2A^-dimensional stochastic differential equations are transformed to 13- 
dimensional deterministic differential equations expressed in terms of means, 
variances and covariances of state variables. The synchronization ratio and 
the firing-time precision for an applied single spike have been discussed as 
functions of Z and p. Our calculations have shown that with increasing p, 
the synchronization is worse because of increased heterogeneous couplings, 
although the average network distance becomes shorter. Results calculated 
by out theory are in good agreement with those by direct simulations. 
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I. INTRODUCTION 



It is well known that a brain forms complex networks with nodes (neurons) and links 
(axons and dendrites). A small patch of cortex may contain thousands of similar neurons, 
each connecting with hundreds or thousands of other neurons in that same patch or in 
other patches through axons and dendrites. The underlying dynamics of individual neu- 
rons is described by Hodgkin-Huxlcy-type nonlinear differential equations (DEs). Many 
theoretical studies have been reported on dynamics of large-scale neuron networks. Exten- 
sive numerical calculations have been made by using various spiking neuron models such as 
Hodgkin-Huxley (HH) [1], FitzHugh-Nagumo (FN) [2,3] and Hindmarsh-Rose (HR) models 
[4]. These theoretical studies have been performed with the use of the two approaches: 
direct simulations and analytical methods such as the Fokker-Planck equation [5], the pop- 
ulation density method [6,7] and the moment method [8-11]. Since the computation time 
of direct simulations is proportional to N"^, simulations for actual network size become diffi- 
cult, where N is the size of a given neuron network. The Fokker-Plank equation method is 
mainly applied to A?^ = cx) network with the mean-field and/or diffusion approximations [12]. 
The population method has been employed for a large-scale integrate-and-fire (IF) neuron 
network [6,7]. The moment method has been applied to FN and HH neuron models [8-11]. 

Most of theoretical studies have assumed that couplings in neuron networks are local 
{Z <^ N) or global {Z — N — 1), and/or regular (p = 0) or random {p — 1), where Z and 
p denotes the average coordination number and the concentration of random couplings, re- 
spectively. In real neuron networks, however, couplings are neither local nor global with the 
degree of randomness locating between the two extremes of regular and random couplings. 
In recent years, much attention has been paid to small-world (SW) networks with the finite 
degree of heterogeneity in couplings, which is characterized by the high clustering and the 
small average distance between nodes [13-16]. The SW property is realized in various kinds 
of biological, social and technological systems such as the electric power grid, the movie-star 
collaborations and the neuronal network of the nematode worm C. elegans [13,14]. Some 
calculations have been reported for neural networks of spiking neuron models as well as of 
phase models [17]- [23]. It has been shown that by introducing the coupling heterogeneity 
into SW networks, the synchronization is increased because the average distance in SW net- 
works is shorter than that in regular networks [17-19] [21-23]. Recently, however, Nishikawa 
et al. [20] have claimed that the synchronization is decreased with including the coupling 
heterogeneity in SW networks. Then it has been controversial whether the synchronization 
in SW networks is better or worse than in regular networks. These studies on SW networks 
have entirely relied on direct simulations, and it is desirable to make a study by using an 
analytical method. 

In previous papers of Refs. [24] and [25] (which are referred to as 1 and 11), the present 
author proposed a semianalytical dynamical mean-field approximation (DMA) theory for a 
study on neuron ensembles (networks) with all-to-all (global) couphngs. In I, DMA was ap- 
plied to an A^-unit FN neuron network, for which 2A^-dimensional stochastic DEs are trans- 
formed to 8-dimensional deterministic DEs expressed by means, variances and covariances 
of state variables. In the following 11, DMA was applied to networks consisting of general 
spiking neurons, each of which is described by M variables. MA^-dimensional stochastic 
DEs are transformed to N^q deterministic DEs where N^q = M{M + 2). The DMA theory 
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was successfully applied to HH neuron network with N^q = 24 in II. Advantages of DMA 
are (i) some qualitative properties of networks are derived without numerical computations, 
and (ii) the computational time of DMA is much shorter than those of the moment method 
[26] and direct simulations. As for the item (ii), for example, the former is thousands times 
faster than the latter for N — 100 HH neuron network with 100 trials [25]. 

The purpose of the present paper is to develop a new approach for SW neural networks 
of FN neurons with general couplings, extending our semianalytical DMA [24] [25]. In I 
and II, interactions among neurons are assumed to be all-to-all (global) couplings. For 
DMA to include local couplings in SW networks, we have taken into account variances and 
covariances which express three kinds of spatial correlations: (1) on-site correlation, (2) the 
correlation for a coupled pair and (3) that for an uncoupled pair without direct couplings. 
Assuming that the heterogeneity is small, we have included its effects in order to discuss the 
synchronization in SW networks. 

The paper is organized as follows. In Sec. II, we have derived DEs, applying the DMA 
to SW networks consisting of FN neurons which are coupled with the average coordination 
number Z. The original 2iV- dimensional stochastic DEs are transformed to 13-dimensional 
deterministic DEs. In Sec. IIIA, we report numerical calculations for regular networks by 
changing Z from local {Z -C N) to global couplings [Z = N — 1). The Z-dependence of the 
firing-time accuracy and the synchronization ratio for an applied single spike is discussed. 
Numerical calculations for SW networks are reported in Sec. IIIB, where the effect of the 
concentration of random couplings is discussed. The final Sec. IV is devoted to conclusion 
and discussion. 

II. SMALL- WORLD NETWORKS OF FN NEURONS 

A. Adopted model and method 

We have assumed that A^-unit FN neurons are distributed on a ring with the average 
coordination number Z and the concentration of random couplings p. Dynamics of a single 
neuron i in a given SW network is described by the non-linear DEs given by 

F[xuit)] -cx,,{t) + lt\t) + lt\t) + 6(t), (1) 
hxu{t)-dx2i{t) + e, {i^lioN) (2) 

I^''\t)^JY.Ci,G{x,,{t)l (3) 
j 

lt\t) ^Ae{t- tin) 0{Un + t^-t). (4) 

kx{t) [x{t)-a] [l-x{t)], k = 0.5, a = 0.1, b = 0.015, d = 0.003 and 
e = [8] [9] [24] : xu and X2i denote the fast (voltage) variable and slow (recovery) variable, 
respectively: G{x) stands for the sigmoid function given by G{x) = 1/(1 + exp[—{x — 9)/a\) 
with threshold 9 and width a: J the coupling strength: Qj the coupling factor given by 



dxii{t) 
dt 

dX2i{t) 

dt 



with 



InEq. (l)-(4),F[a;(t)] 
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Cij — Cji = 1 for a coupled (i, j) pair and zero otherwise, self-coupling terms being excluded 
(cii = 0). By changing Z value, our model given by Eqs. (l)-(4) covers from local couplings 
{Z <^ N) to global couplings {Z = N — 1). Wc have studied the response of neuron 
networks to an external, single spike input given by l'f\t) with magnitude A and spike 
width tyj applied at the input time Q{x) being the Heaviside function. Added white 
noises ^^(i) are given by 

< m > = 0, (5) 

<mij{t') > =/5'%5(t-t'), (6) 
where the average of < C/(z, > for an arbitrary function of U{z, t) is given by 

<U{z,t) >= J ... J dzU{z,t) Pr{z), (7) 

Pr{z) denoting a probability distribution function for 2A^-dimensional random variables 

z = {{x^i}). 

An SW network is made after the Watts- Strogatz model [13]. Starting from the regular 
coupling for which = Cojj, Nch couplings among NZ/2 couplings are randomly modified 
such that CQij = is changed to = 1 or verse versa. The concentration of random 
couplings is given by 

which is and 1 for completely regular and random couplings, respectively. We shall take 
into account the effect of the heterogeneity given by 



= — c 



Z Z 



coij), (9) 



assuming it is small. 

After I, we will obtain equations of motions for means, variances and covariances of state 
variables. Variables spatially averaged over the ensemble are defined by 



and their means by 



^-(t) = ^ E = 1' 2 (10) 



^,,{t)^{{x,{t)))^, (11) 



where the bracket < ■ >c denotes the average over the coupling configuration. As for 
variances and covariances of state variables, we consider three kinds of spatial correlations: 
(1) on-site correlation (7), (2) the correlation for a coupled pair (C) and (3) that for a pair 
without direct couplings [r]): 



' 1k,x, for i = j 

C«,A, for i ^ j, Cij = 1 (12) 
, Vk,\: for i 7^ J, Cij = 



4 



where k,,X—1, 2 and 



(13) 



In Eq. (12), 7,,^, Ck,x 



and 77k,a a-re defined by 




c 



(15) 



(14) 




Cij) {Sx^i{t) Sxxj{t)) 



(16) 



For a later purpose, we define also the spatially-averaged correlation given by 




(17) 



{{5X^{t) 5X,m^, 



(18) 



where 6X^(t) = X^(t) — ji^it). It is noted that 7^;^^; Ck,a ^k,a and p^.A are not independent, 
obeying the sum rule given by 



with Cii — Q. 

In calculating means, variances and covariances given by Eqs. (11) and (14)-(17), we 
have assumed that (1) the noise intensity is weak, (2) the distribution of state variables 
takes the Gaussian form, and (3) the coupling heterogeneity of 5cij/Z is small. By using 
the first assumption, we expand DEs given by Eqs. (l)-(4) in a power series of fiuctuations 
around means. The second assumption may be justified by some numerical calculations for 
FN [9] [27] and HH neuron models [28] [29]. Based on the third assumption, the effect of 
coupling fluctuations has been taken into account up to the order of 0{{6cij/ ZY). 

Before closing Sec. IIA, we briefly summarize the introduced variables and their mean- 
ings as follows: A^, the number of neurons: Z, the average coordination number: p, the 
concentration of random couplings: J, the coupling strength: Qj, the coupling factor be- 
tween neurons i and j: X^, the spatially average of the fast {k — 1) and slow {k — 2) 
variables; a mean value of X^] ^k,X) Ck,x, and ry^.A; the correlations of on-site, a coupled 
pair and an imcoupled pair, respectively. Readers who are not interested in mathematical 
details, may skip to Sec. IIC where a summary of our method is presented. 



Np^,x = 7«,A + ZCn,x + {N-Z- 1)?7„,A- 



(19) 



In order to derive Eqs. (14)-(19), we have employed the decomposition: 



1 = Sij + {1- Sij)[cij + (1 - Cij)], 



(20) 
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B. Equations of motions 

After some manipulations, we get the following DEs (the argument t being suppressed; 
for details, see appendix A): 



~dt 



fo + /27i,i - c/i2 + JZ{go + gi(f)i) + lext, (21) 
bui - dfX2 + e, (22) 



^ = 2(a7i,i - C7i,2) + 2 JZ((7iCi,i + ^o^i) + (23) 
dt 

^7l,2 

dt 

dpi,i , , f2JZgi\^_ , , ^„ ^ , 

dt 
dp2,2 



2(671,2-^^72,2), (24) 
671,1 + (a - rf)7i,2 - C72,2 + JZ{gi(i^2 + 5002), (25) 



2(api,i - cpi,2) + [-^) [71,1 + ZRCi,, + (N-ZR- 1)771,1] + ^, (26) 

-2(6pi,2-dp2,2), (27) 

^ = 6pi,i + (a - d)pi,2 - cp2,2 + (^) [71,2 + ^i?Ci,2 + (AT - Zi? - 1)771,2], (28) 

^ = 2(aCi,i - cCi,2) + 2 Jyi[7i,i + ^CCi.i + {ZR - ZC - l)?7i,i], (29) 

^ = 2(6Ci,2-ciC2,2), (30) 

^ = 6Ci,i + (a - d)Ci,2 - cC2,2 + J5'i[7i,2 + ^CCi,2 + {ZR - ZC - 1)771,2], (31) 
= 2(a77i,i - C77i,2) 

+ [ n-^Z- 1 ) ^^'f^^^ " " ^^^^'^ + " + ^^^"''''^ " ^"'^'^ ' ^^^^ 

= 2(6771,2 - (^772,2), (33) 

= 67/1,1 + (a - d)77i,2 - C7/2,2 

+ ( ^/| J {5i[(^i? -ZC- l)Ci,2 + (A^ - 2Zit: + ZC)77i,2] - 5002} , (34) 

^^■^ ^ acpi - c(l)2 + JZgodRp, (35) 

(36) 



— - = b(pi - d(p2, 

dt 



with 



Mt) = EE (^^^^c,,)) , = 1, 2 (37) 
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C — X! X! X] coijCojfcCoife, (38) 

i i h 

^ = ^oij cojk, (39) 

i j k 
\ i j k / ^ 

where a = /i + 3/371,1, fe = {l/i\)F^^\ gi = {l/i\)G^^\ C corresponds to the clustering 
coefficient introduced in SW networks [13,14], R expresses the couphng connectivity, and 
5Rp is its fluctuation part, related discussions being given in Sec. IV. 



C. Summary of our method 

Equations of motions for ii^{t), 'y^A^), Ck,a(^), VK,\{t) and p^A^) are given by Eqs. (21)- 
(40). In Eqs. (35) and (36), (pnit) {k — 1, 2) are new correlation functions which appear 
in the process of calculating equations of motion of 7k,a et al. The factors C, R and SRp 
deflned by Eqs. (38)- (40) generally depend on the geometry of a given neuron network. For 
a regular ring with even Z, we get R—1 and C given by 



C 



0, for Z < 2 

3/4-3/2Z, for4<Z<2A^/3 

3/4 - 3/2Z + 9/4 - {3N - 9/2) /Z (41) 

+ {N^ -3N + 2)/Z^, for 2iV/3 <Z<N-1 

(1-1/Z). iorZ = N-l 



Figure 1 shows C as a function of Z/N for N = 100, 200, 500 and 1000. Wc note that 
C ~ 0.75 for 0.1 < Z/N < 0.7 and that C ^ {1 ~ 1/Z) as Z/N ^ (1 - 1/N). In the 
case of global couplings {Z = N — 1), however, we get C = (1 — l/Z) independent of the 
geometry. 5Rp deflned by Eq. (40), which expresses fluctuations in heterogeneous couplings, 
is increased with increasing the concentration of random couplings, p [Fig. 6(a)]. Among 
the 12 correlations such as 7k,a al. given by Eqs. (14)-(17), 9 correlations are independent 
because of the sum rule given by Eq. (20). In this study, we have chosen nine correlations 
of 7k,a, Ck,a and p^.A as independent variables. Then the original 2iV-dimensional stochastic 
DEs given by Eqs.(l) and (2) have been transformed to 13-dimensional deterministic DEs. 

It is worthwhile to explain the relation between the present theory and I where the 
original 2A^-dimensional stochastic DEs for regular, global couplings are transformed to 8- 
dimensional deterministic DEs. In the present study for the general coupling, we have to 
take into account (k,x and rj^^x, in order to discriminate correlations between a coupled pair 
and an uncoupled pair. However, in the limit oi Z — N — 1 ior regular, global couplings for 
which R — 1 and ZC — Z — 1, rj^^x are not necessary because there are no uncoupled pairs: 
prefactors of {ZR — ZC — 1) for rji^ x in Eqs. (32) and (34) vanish with 0« = 0. Then the 
number of required DEs is reduced from 13 to 8. Equations (21)- (28) for 7k;,a and p^^A 
agree with Eqs. (20)-(27) in I [30]. 
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D. Firing-time accuracy and synchronization 

Firing-time accuracy 

When we solve DEs given by Eqs. (21)- (36), we may obtain various quantities relevant 
to firings in neuron networks. The firing time of a given neuron i is defined as the time when 
the variable xii{t) crosses the threshold 9 from below: 

toe = {t\xuit) = e;x,{t)>0}. (42) 

It has been shown that the distribution of firing times of tot is given by [24] 

' ^ (43) 



with 




V7i,i(^/) 

for 7i,i(t/) 




= ^ , (44) 

where $ expresses the normal distribution function, the average firing time tf is implicitly 
defined by [Xiitf) = 9, fii = fii{tf) and the dot denotes the time derivative. 

Similarly, the firing time of an averaged variable Xi{t) is defined as the time when the 
variable Xi{t) crosses the threshold 9 from below: 

tog^{t\X^{t)^9;Xi{t)>0}. (45) 

The distribution of firing times of tgg is given by [24] 

Z.it) - * i 1 0(Ai), (46) 



with 




5{t-tf). forpi,i(t/) 



Kg^Jf^. (47) 



Synchronization ratio 

We discuss the synchronization in neuron networks, considering the quantity given by 

^^W = 4EE<M^)-^.W]'>> (48) 

* j 

= 2(71,1 -pi,i), (49) 

which vanishes in the completely synchronous state. Prom a comparison of Eqs. (23)- (25) 
with Eqs. (26)- (28), we note that 
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P«,A=^, forJ^O (50) 

Then, Rs{t) given by Eq. (49) becomes Rs{t) = {1 — 1/N)'~fi i{t) = Rso{t) in the asynchronous 
state, while Rs{t) = in the completely synchronous state. We define the synchronization 
ratio at the firing time t/ by [24] 

Sf^Sitf), (51) 

with 

which is and 1 for completely asynchronous and synchronous states, respectively. The syn- 
chronization ratio shows much variety depending on model parameters such as the coupling 
strength (J), the noise intensity (/3), the size of cluster (N), the coordination number (Z), 
and the random concentration (p), as will be discussed in Sec. III. 



III. CALCULATED RESULTS 
A. Regular couplings 

We have adopted same parameters of ^ = 0.5, a — 0.5, — 10, A — 0.10, tin = 100 and 
Tu, — 10 as in I [24]. DMA calculations have been made by solving Eqs. (21)-(36) with the 
use of the fourth-order Runge-Kutta method with the time step of 0.01. We have performed 
direct simulations by using also the fourth-order Rungc-Kutta method with the time step 
of 0.01. Results of direct simulations are averages of 1000 trials for Z < 20 (or < 20) and 
those of 100 trials otherwise noticed. All quantities are dimensionless. 

First we discuss the case of regular coupfings {p — 0), by changing the average coordi- 
nation number Z from local {Z <^ N) to global couplings (Z = N — 1). The plots in Figs. 
2(a)-2(c) show firings in an = 100 neuron network with regular couplings for Z = 10, 
50 and 99 with P = 0.01 and J = 0.002 when a single external spike given by Eq. (4) is 
applied. Figures 2(a)-2(c) show that as increasing Z, scattering of firing times is reduced, 
which suggests an increase in the firing accuracy and the synchronization. These are results 
of direct simulations with single trials. They are more clearly discussed with calculations 
using the DMA. Fig. 2(d)-2(f) show time courses of S{t) calculated in the DMA for Z = 10, 
50 and 99, whose magnitudes are increased as increasing Z; note differences of the ordinate 
scales in Figs. 2(d)-2(f). The synchronization ratio at firing times, Sf, is 0.0019, 0.0113 
and 0.0295 for Z — 10, 50 and 99, respectively, which shows an increased synchrony with 
increasing Z. 

We will discuss some details of the DMA calculation in Figs. 3(a)-3(d) which show time 
courses of /xi, 71^1, ^1^1, and pi^i, respectively, for regular couplings {p = 0) with P = 0.01, 
J = 0.002, N = 100 and Z = 10. Results of DMA expressed by solid curves are in good 
agreement with those of direct simulations depicted by dashed curves. Time courses of 
/^i; 7i,i and pi^i shown in Fig. 3(a), 3(b) and 3(d) for local coupfings (Z=10) are not so 
different from those for global couplings having been reported in Fig. 1 of I, except for their 
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magnitudes. For example, DMA calculations for the local coupling with Z = 10 in the case 
of (3 = 0.01, J = 0.002 and N = 100 show that magnitudes of 7ij, and pi i at the firing 
time of t = 104.44 are 0.271 xlO'^, 0.475 xlO"'' and 0.320 xlO'^ respectively. In contrast, 
for the global coupling with Z — 99, magnitudes of 71^1, and pi^i at the firing time of 
t = 103.88 are 0.235 xlO-^ 0.693 xlO"^ and 0.921 x 'lO"^, respectively. 

Figure 4(a) shows the Z dependence of 71,1, and at the firing time with J = 
0.002, f3 = 0.01 and N = 100; filled and open marks express results of DMA and direct 
simulations, respectively. Results of 71^1 and pi^i of DMA are indistinguishable from those 
of direct simulations. With increasing Z, both and pi^i are increased, while 71^1 is 
slightly decreased, as mentioned above. The Z dependence of the firing time tf is plotted 
in Fig. 4(b), which shows the faster response for larger Z. This is due to the fact that by 
an increased Z, fii is increased more rapidly to cross the threshold level of 6. Then fii at 
t = tf is increased with increasing Z, as the chain curve in Fig. 4(c) shows. Figure 4(c) 
shows that with increasing Z, the firing-time accuracy of Stot is improved while that of Stgg is 
independent of Z. The Z dependence of the synchronization is plotted in Fig. 4(d) showing 
Sf to be linearly increased for a small Z. This clearly explains the larger synchrony Sf for 
larger Z, having been shown in Figs. 2(a)-2(f). 

B. SW couplings 

Next we discuss the case of SW couplings, by changing the concentration of random 
couphngs p. The plots in Figs. 5(a)-5(c) show firings in SW networks for p — 0.0, 0.1 
and 1.0, respectively, with (3 = 0.005 J = 0.02, = 100 and Z = 10 calculated by direct 
simulations with single trials, when a single external spike given by Eq. (4) is applied. In this 
subsection, we have adopted a smaller /? and a larger J than in Sec. IIIA to get more evident 
effects of p. Figures 5(a)-5(c) show that as increasing p, scattering of firing times is gradually 
increased, which suggests a decrease in the firing accuracy and the synchronization. These 
results are more clearly seen in calculations with the use of DMA. Fig. 5(d)-5(f) show time 
courses of S{t) for p = 0, 0.1 and 1.0, calculated in the DMA. The synchronization ratio at 
firing times Sf is 0.0256, 0.0224, and 0.0114, for p = 0, 0.1 and 1.0, respectively. Although 
Sf for p = 0.1 is nearly equal to that for p = 0.0, the time course of S{t) for p = 0.1 is rather 
different from that for p = 0.0. 

This decrease in 5"/ with increasing p mainly arises from an increased 5Rp, as shown in 
Fig. 6(a) where the p dependence of SRp is plotted for Z = 10, 20 and 50 of a given ring 
with N = 100. With increasing p, 6Rp is linearly increased as SR oc p/Z for a small p. 
Figure 6(b) will be explained in Sec. IV. 

Figure 7(a) shows the p dependence of 71,1, and pi^i at the firing time with J = 0.02, 
P — 0.005, N — 100 and Z — 10; filled and open marks express results of DMA and direct 
simulations, respectively. At p = 0.0, 71^1, and pi^i are 0.671 x 10~^, 0.131 x 10~^ and 
0.239 X 10""^, respectively. In contrast, at p = 1.0, they are 0.109 x 10"^ 0.144 x 10"^ 
and 0.232 x 10~^, respectively. With increasing p, 71^1 is increased, while pi,i and are 
almost constant. The difference between the p dependences of 71,1, and arises from 
the fact that dji^i/dt in Eq. (23) has a contribution from 0i while dpi^i/dt and d(i^i/dt in 
Eqs. (26) and (29) have no direct contributions from it. Figure 7(b) shows that the firing 
time oftf — 103.88 is independent of p, which is in accordance with a constant fii shown in 
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Fig. 7(c). Figure 7(c) shows that with increasing p, the firing-time accuracy of Stoe becomes 
worse because of an increased 71^1 while that of 6tog is independent of p. The p dependence 
of Sf is depicted in Fig. 7(d), which shows that the synchrony is decreased with increasing p. 
This clearly explains results of smaller Sf for larger p, having been shown in Figs. 5(a)-5(f). 



IV. CONCLUSION AND DISCUSSION 

Generalizing a phenomenological analysis adopted in I [24] based on calculated results of 
DMA, we have tried to get an analytical expression for Sf. From calculated results discussed 
in the previous section, we expand 71,1 and in a series of JZ and p: 

ji,i^^o[i--aiJZ{l-a2p) + --], (53) 
p,,,^^(l + hJZ + ..), (54) 

where 70 oc and ai, 02 and bi are positive coefficients. We have obtained an expression 
for 71^1 given by Eqs. (53), because the effect of p should vanish for J = or Z = 0. 
Substituting Eqs. (53) and (54) to Eq. (52), we get 

S,= [^i^^^^li^)jZ^... (55) 

The expression for 5"/ given by Eq.(55) well explains the behavior shown in Figs. 4(d) and 
7(d). Dependences of the quantities on N, Z, J and (3 for local couplings are the same as 
those for all-to-all couplings having discussed in I. Typical examples of dependence of 
various quantities are shown in Figs. 9(a)-9(d). Figures 9(a) and 9(b) show that pi^i oc 
while 7i,i, and tf are independent of N, which yields Stoe oc N~^/'^ and Stag oc N^, as 
shown in Fig. 9(c). Figure 9(d) shows that Sf ex both for local and global couphngs, 
expressing that the synchronization is more easily realized in smaller networks than in larger 
ones. 

In an early stage of this study, we obtained DEs given by Eqs. (21)-(34) with (pi = 4>2 = 0, 
but with C and R which are replaced by Cp and Rp, respectively, given by [for detail see 
after Eq. (A22) in appendix A] 

Cp ^ CijCjkCik^ , (56) 

= ( E E E cjk) ■ (57) 

\ i j k If. 

In this formulation, the effect of the couplings heterogeneity is included in the p-dependent 
clustering coefficient Cp and coupling connectivity Rp. The clustering coefficient Cp denotes 
an averaged fraction for given three nodes to be mutually coupled [13,14]. The p dependence 
of Cp is depicted in Fig, 6(b) which shows that with increasing p, Cp is decreased and 
approaches Cp — Z/N at p — 1. In contrast, the coupling connectivity Rp expresses an 
averaged fraction for given two nodes, which are not necessarily coupled, to have a common 
neighboring node. Rp in Eq. (57) may be rewritten as 
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^P = ^ E^'^W^^^'' (58) 

where the overhne denotes the average over P{K) expressing the probabihty for a given 
neuron to have K couphngs [31]. It is easy to see that Rp is given by i?p = 1 + 5Rp [Eqs. 
(60) and (61)], the p dependence of 5Rp being plotted in Fig. 6(a). Unfortunately results 
calculated with the use of Cp and Rp for finite p were not in good agreement with those 
of direct simulations because effects of coupling heterogeneity are not properly taken into 
account in such DEs. 

After several tries, we have obtained DEs having been given by Eqs. (21)-(36). C, R 
and 5Cp given by Eqs. (38)-(40) may be expressed in terms of Cp and Rp as [31] 

C = Co, (59) 
i? = i?o = 1, (60) 

5Rp = Rp-Ro = ^{K-K)\ (61) 

with Z ^K. Figure 8 shows P{K) for p = 0.0, 0.1, 0.2 and 1.0 with N = 100 and Z = 10. 
In the limit of p = 0, P{K) {— 5k,z) is the delta function. With increasing p, P{K) has 
the distribution centered dX K = Z . In the limit of p = 1, P{K) approaches the Poisson 
distribution [16]. Figure 6(a) shows that with increasing p, 6Rp is increased, while Cp is 
decreased as shown in Fig. 6(b). An increased 5Rp yields an increase in 71 1, by which Sf is 
decreased and Stoe is increased. It should be noted that effects of heterogeneous couplings 
are taken into account by 5Rp through the correlation functions 0i and 02 in Eqs. (35) and 
(36), which play important roles in dynamics of SW networks. 

To summarize, we have developed a semianalytical theory for SW networks of spiking 
FN neurons, including three kinds spatial correlations: correlations of on-site, a coupled pair 
and an uncoupled pair. By changing Z and p, we have performed model calculations of the 
response of the network to an external single spike. It has been shown that 

(1) when Z is increased, the synchronization ratio Sf and the firing-time accuracy Stoe are 
improved [Figs. 4(c) and 4(d)], which arises from a decrease in 71^1 and an increase in pi^i, 
and 

(2) when p is increased, both Sf and 6toe become worse [Figs. 7(c) and 7(d)] due to an 
increase in 71 1 induced by fiuctuations in the coupling heterogeneity. 

The item (1) is easily understood. The result for Sf in the item (2) is consistent with 
that of Ref. [20]. It, however, contradicts some calculations [17-19] [21-23], which have 
claimed that the synchronization in SW networks is better than that in regular networks, 
since communication between neurons is more efficient because of the shorter characteristic 
path length L (as for the p-dependence of L, see Fig. 2 of Ref. [13]). Our semianalytical 
theory with the use of the DMA, which is valid for weak noise (/9 ^ 1) and small coupling 
heterogeneity {6Rp <^ 1), has shown that the synchrony of SW networks depends on R, 
C and SRp given by Eqs. (38)- (40), but it is not affected by the average path length, 
L. In particular, 6Rp, cpi and 02 have been shown to play crucial roles in dynamics of 
SW neural networks. Although the item (2) discussed above rehes on the definition of the 
synchronization ratio of S(t) given by Eq. (52), this conclusion is not changed even if we 
adopt an alternative measure for the synchrony. For example, when we employ Rg given by 
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Eq. (49), Rs is increased with increasing p because of an increased 71,1, which again signifies 
the worse synchronization in SW networks. The semianalytical theory developed in this 
paper can be apphed not only to SW neural networks but also to a wide class of complex 
SW networks. When we apply our theory to a general SW network in which dynamics of each 
node is described by M-dimensional stochastic DEs, we get A^eg-dimensional deterministic 
DEs where Neq = M{3M + 7)/2. For example, N^q = 5 for Langevin model (M = 1), 
Neq = 13 for FN model (M = 2), and N^q = 38 for HH model (M = 4). The items (1) and 
(2) [and also Eq. (55)] which have been derived for FN neuron model, are expected to hold 
for any SW network. 

The present approach shares in its advantages with the original DMA previously pro- 
posed in I: (i) some results may be derived without numerical calculations because of its 
semianalytical nature, and (ii) a computational time for a large-scale system by DMA is 
much shorter than that by direct simulations. By extending the ring geometry adopted in 
this paper, we may discuss the response of more realistic synfire-chain-type SW networks 
[24] [32] . In the present paper, we have neglected the transmission time delay. Because the 
average path length L becomes shorter by the appearance of shortcuts [13-16], the response 
speed is expected to be improved in SW networks with time delays. Recently, we have 
successfully applied the DMA to stochastic ensembles with time-delayed regular couplings 
[33,34]. By using our approach, we may discuss dynamics of general SW networks with time 
delays within the framework of the DMA. In the so-called scale-free (SF) network such as 
the world-wide web and the network of citations of scientific papers, the link connectivity 
P{K) for a node to interact to K other nodes follows a power-law distribution P{K) ~ K'"' 
with the index 7 (~ 2.1 to 4) [35], in contrast to an exponential distribution for a large 
K in our SW networks (Fig. 8). This SF distribution probability originates from the two 
factors, the growth of nodes and their preferential attachment [35]. Quite recently it has 
been reported that the functional connectivity P{K) versus the distance K in human brain 
is given by a SF distribution: P{K) ~ K^^ [36]. It is interesting to apply our semianalytical 
approach to such SF networks. These subjects raised above are left as our future study. 
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Substituting Eqs. (9) and (13) to Eqs. (l)-(4), we get DEs for 5xu and 5x2i of a neuron 
i, given by (argument t is suppressed) 
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APPENDIX A: DERIVATION OF EQS. (21)-(36) 



d5x 



li 



fiSxu + f2{Sxli - 7i,i) + fsSxfi - c6x2i + SIl 



(Al) 



d.t 



d6x2j 



b5xij — d5x2j, 



(A2) 



dt 



with 
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5I^''\t) = J Y^[gi(t)coij5xij(t) + go{t)5cij + gi(t)5cij5xij(t) 

.7 



(A3) 



where = (1/^!)FW and = (l/f!)^^. DEs for the correlations are given by 



dt \ Ar ^ 



dt 



EE' 



dt 



+ 



^ d6xK 



dt 



d6x 



+ 



KJ 



dt 



Sxm 



( /v2 EE 



' d5x 



dt 



'd_5x^\ 



Xi 



(A4) 
(A5) 
(A6) 



With the use of Eqs. ( Al)-( A3), wc may calculate DEs given by Eqs. (21)-(34). For example, 
terms including Slj;'^^ in d'ji^i/dt, d(i^i/dt and dpi^i/dt become 



N 

2, J 



i j 



\ i j 



+ ^EE5o((5a;H5Qj))^, 

i j 

2J 

= ]^ E E E 9iCoijCojk {{SxiiSxik))^ 



i j k 



-i, EE (^-1.^4 



-(c) 



+ E E E ^oCoii ((^Xii^Cjfc))^ , 

1' J k 

= 2Jyi[7i,i + ZCCi,i + (ZR -ZC- 1)771,1], 
2J 

= ■^Y.J2J29iCojk{{Sxii6xik))c 

i j k 

2J 

+ l^EEEfo(('5a;ii(59fc))^, 

i j k 

= ^^[71,1 + ^^Cl,l + (AT - Zi? - 1)771,1], 



where 0^ = 1, 2) are new correlation functions defined by 



K=l, 2 

c 

In evaluating Eqs. (A7)-(A12), we have employed the relations 



given by 



1 = 



1 



NZ^, 



EE'^oij, 



2 J 



1 



NZ^ 



E E E '^oij cojk, 

i j k 



(AT) 
(A8) 

(A9) 
(AlO) 

(All) 
(A12) 

(A13) 



(A14) 
(A15) 
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C - XI 51 XI '^oij cojk coik, (A16) 

i j k 

= ]^ (E E E Sc., 5c,^ , (A17) 

and the mean-field approximation given by 

{{Sx^i5 xxj))^ = 7«,A Sij + (1 - 6ij) [Cn,x Sij Cij + r]^^x 5,j (1 - c^)], (A18) 

= 7/c,A ^ii + Ck,A + 77«,A (1 - % - Cij), (A19) 

{{5x^^i5cjk))^ = 0«Cjfc (5y + (A20) 

with the Gaussian decouphng approximations [24]. In Eqs. (A18) and (A19), 7k,a, Ck,a 
and ?7k,a denote the correlations of on-site, a coupled pair and an uncoupled pair, which are 
defined by Eqs. (12)-(14). The approximations given by Eqs. (A18)-(A20) are consistent 
with the definition of jf^x, (f^x and rj^x given by Eqs. (14)-(16), and those of 4>^ given by 
Eq. (37). 

The equations of motion of 0^ are similarly calculated with the use of the relation given 

by 

which yield Eqs. (35) and (36). 

We have taken into account terms up to orders of 0{{6xY), 0{{6c/Z)'^) and 0{6x 6c/Z) 
in Eqs. (21)-(36), and up to the order of 0{{6x)'^) in the term including a {= fi + 3/371,1) 
which plays an important role in stabilizing DEs [24]. 

On the contrary, when we adopt an expression given by 

Sll'\t) = J J2[g^it)c,,dx,jit) + •], (A22) 
j 

instead of Eq. (A3), DEs given by Eqs. (A7), (A9) and (All) become 

(|E('5^h^^^)) ^^EE5i(c..(('5^A.))eX, (A23) 

\ i / c i 3 

- 2JZ^iCi,i, (A24) 

(2 \ U 

J^Y^Y. % {5xu5lf)) - ^ E E E ^1 (cyCjfc {{SxiiSxik)) X ' (A25) 



2J(7i[7i,i + ZC/1,1 + {ZRp - ZCp - 1)7^1,1], (A26) 



W2YY (SxuSif)) ^ 1^ E E E ^1 {cjk {{SxuSx,,))X , (A27) 



^^[71,1 + ZR^Ci,! + {N- ZR, - 1)771,1], (A28) 



where decoupling approximations such as 
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{cij {5xii 5xij))^ ~ (^Cij {{5xii Sxij))^^^ , (A29) 

and Eq. (A19) are employed. Cp and Rp in Eqs. (A26) and (A28) are given by Eqs. (56) 
and (57). Note that Cij in Eqs. (A23), (A25) and (A27) depends on the configuration of 
couphngs while coij in Eqs. (A7), (A9) and (All) does not. Then we got equations of 
motions given by Eqs. (21)-(34) with 0i = 02 = but with C and R which are, respectively, 
replaced by p dependent Cp and Rp given by Eqs. (56) and (57). As mentioned in Sec. IV, 
results calculated with the use of such DEs arc not in good agreement with those obtained 
by direct simulations because effects of coupling fiuctuations are not properly included in the 
formulation mentioned above. It is indispensable to take into account effects of the coupling 
heterogeneity expressed by SRp through the correlation functions 0i and 02, as given by 
Eqs. (35) and (36). 



16 



REFERENCES 



[1] A. L. Hodgkin and A. F. Huxley, J. Physiol. 117, 500 (1952). 
[2] R. FitzHugh, Biophys. J. 1, 445 (1961). 

[3] J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962). 
[4] J. L. Hindmarsh and R. M. Rose, Nature 296, 162 (1982). 

[5] H. Risken, The Fokker- Planck Equation: Methods of Solution and Applications, Springer 

Series in Synergetics, Vol. 18 (Springer Verlag, 1992). 
[6] B. W. Knight, J. Gen. Physiol. 59, 734 (1972); B. W. Knight, Neural Comput. 12, 473 

(2000); B. W. Knight, A. Omurtag, and L. Sirovich, Neural Comput. 12, 1045 (2000); 

A. Omurtag, B. W. Knight, and L. Sirovich J. Comput. Neurosci. 8, 51 (2000). 
[7] D. Q. Nykamp and D. Tranchina, J. Comput. Neurosci. 8, 19 (2000); E. Haskell, D. Q. 

Nykamp, and D. Tranchina, Network 12, 141 (2000). 
[8] R. Rodriguez and H. C. Tuckwell, Phys. Rev. E 54, 5585 (1996). 
[9] H. C. Tuckwell and R. Rodriguez, J. Comput. Neurosci. 5, 91 (1998). 
[10] R. Rodriguez and H. C. Tuckwell, BioSytems 48, 187 (1998). 

[11] R. Rodriguez and H. C. Tuckwell, Mathematical and Computer Modeling 31, 175 
(2000). 

[12] N. Fourcaud and N. Brunei, Neural Comut. 14, 2057 (2002). 

[13] D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998). 

[14] S. H. Strogatz, Nature 410, 268 (2001). 

[15] M. E. J. Newman, J. Stat. Phys. 101, 819 (2000). 

[16] R. Albert and A. Barabasi, Rev. Mod. Phys. 74, 47 (2002). 

[17] L. F. Lago-Fernandez, R. Huerta, F. Corbacho, and J. A. Sigiienza, Phys. Rev. Lett. 
84, 2758 (2000). 

[18] M. Barahona and L. M. Pecora, Phys. Rev. Lett. 89, 54101 (2002). 
[19] M. Bucolo, L. Fortuna, M. La Rosa, Chaos, Solitons and Fractals 14, 1059 (2002). 
[20] T. Nishikawa, A. E. Motter, Y. Lai, and F. C. Hoppensteadt, Phys. Rev. Lett. 91, 14101 
(2003). 

[21] O. Kwon and H. Moon, Phys. Lett. A 298, 319 (2002). 

[22] H. Hong, M. Y. Choi, B. J. Kim, Phys. Rev. E 65, 26139 (2002). 

[23] H. Hong, M. Y. Choi, B. J. Kim, Phys. Rev. E 65, 47104 (2002). 

[24] H. Hasegawa, Phys. Rev E 67, 041903 (2003). 

[25] H. Hasegawa, Phys. Rev E 68, 041909 (2003). 

[26] When the DMA is applied to A^-unit neuron networks with global, regular couplings, in 
which dynamics of each neuron is described by M variables, MA^- dimensional stochastic 
DEs are replaced by A'eq- dimensional deterministic DEs where Neq = M[M + 2) (Ref. 
[25]). In contrast, when the moment method (Ref. [8]) is employed, we get N^q = 
{1/2)MN{MN + 3). In the case of = 1000, for example, we get N^g = 8 (2 003 000) 
for M = 2, and 24 (8 006 000) for M = 4 in the DMA (the moment method). 

[27] S. Tanabc and K. Pakdaman, Phys. Rev. E 63, 31911 (2001). 

[28] S. Tanabc, S. Sato, and K. Pakdaman, Phys. Rev. E 60, 7235 (1999). 

[29] S. Tanabe and K. Pakdaman, Biological Cybernetics 85, 269 (2001). 

[30] The normalization factor in front of the coupling term is in I while it is unity 
in this paper; results of the former are obtainable from the latter by a replacement of 
J^w/N. 



17 



[31] We may show that Z =< (l/N) EijCij >c= T.kKP{K) and 

Rp =< {l/NZ^)Y.ijkCijCjk >c= Z-^J2kK'^P{K), by defining the probability: 
P{K) =< (l/N) 5iK - Q,) >e. 

This implies that = Rp~Ro = Z''^Y.k K'^[P{K)-Po{K)] = Z-^ K^P{K)-1 = 
Z'^ T.K{K-ZfP{K) which leads to Eq. (61), Pq{K) (= 5k,z) denoting the probability 
for j9 = 0. 

[32] M. Abclcs, H. Bergman, E. Margaht, and E. Vaadia, J. Neurophys. 70, 1629 (1993). 

[33] H. Hasegawa, Phys. Rev. E 70, 021911 (2004). 

[34] H. Hasegawa, Phys. Rev. E 70, 021912 (2004). 

[35] A. Barabasi and R. Albert, Nature 286, 509 (1999). 

[36] D. R. Chialvo, E-print: cond-mat/0402538. 



18 



FIGURES 



FIG. 1. The clustering coefficient C for a ring with regular couplings (jp = 0) as a function of 
Z/N for = 100, 200, 500 and 1000. 

FIG. 2. (color online). The plots showing firings in a regular neuron network for (a) Z = 10, 
(b) 50 and (c) 99 calculated by direct simulations (single trials), and time courses of S{t) for (d) 
Z = 10, (e) 50 and (f) 99 calculated by DMA (solid curves) and direct simulations (dashed curves) 
{(3 = 0.01, J = 0.002, N = 100 and p = 0.0). Arrows in (d)-(f) denote firing times. 

FIG. 3. (color online). Time courses of (a) fii, (b) 71,1, (c) and (d) for j3 = 0.01, 
J = 0.002, N = 100, Z = 10 and p = 0, solid and dashed curves denoting results of DMA and 
direct simulations, respectively. At the bottom of (a), an input signal is plotted. 

FIG. 4. (color online). The Z dependence of (a) the correlations of 71^1 (circles), ^i^i (triangles) 
and pi_i (squares) at the firing time, (b) the firing times tf, (c) the firing-time accuracy of 6toe 
(circles) and 5tog (squares) and fii (triangles), and (d) the synchronization ratio at the firing time, 
Sf, for j3 = 0.01, J = 0.002 and N = 100: filled and open marks denote results of DMA and direct 
simulations, respectively. Results of and fii are only for DMA. 

FIG. 5. (color online). The plots showing firings in a small-world neuron network for (a) 

p = 0.0, (b) 0.1 and (c) 1.0 calculated by direct simulations (single trials), and time courses of 
S{t) for (d) p = 0.0, (e) 0.1 and (f) 1.0 calculated by DMA (solid curves) and direct simulations 
(dashed curves) (/3 = 0.005, J = 0.02 and N = 100). Arrows in (d)-(f) denote firing times. 

FIG. 6. Thep dependence of (a) the factor dRp and (b) the clustering coefficient Cp, for Z = 10, 
20 and 50 with N = 100. 

FIG. 7. (color online). The p dependence of (a) the correlations of 71,1 (circles), (triangles) 

and pi^i (squares) at the firing time, (b) the firing times tf, (c) the firing-time accuracy of Stot 
(circles) and 6tog (squares), and fii (triangles), and (d) the synchronization ratio at the firing time, 
Sf, for P = 0.005, J = 0.02, N = 100 and Z = 10: filled and open marks denote resuhs of DMA 
and direct simulations, respectively. Results of and /ii are only for DMA. 

FIG. 8. The probability P{K) for a given neuron to have K couplings for p = 0.0, 0.1, 0.2 and 
1.0 with N = 100 and Z = 10 in a SW ring . 

FIG. 9. (color online). The N dependence of (a) the correlations of 71,1 (circles), Ci^i (triangles) 

and pi^i (squares) at the firing time, (b) the firing times tf, (c) the firing-time accuracy of St^i 
(circles) and 5tog (squares), and (d) the synchronization ratio at the firing time, Sf {[3 = 0.01, 
J = 0.002, A'^ = 100 and Z = 10): filled and open marks denote results of DMA and direct 
simulations, respectively. In (d) results for global couplings [Z = N — 1) are also shown. 
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This figure "figl-4.gif" is available in "gif" format from: 
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